#include"myownlib.h"

#ifdef MY_USE_QT


std::ostream & operator << (std::ostream & out,QString qstring){
    out<<qstring.toStdString();

    return out;
}


QString getchoosefile(QWidget * parent,string files,QString dir){
    return  QFileDialog::getOpenFileName(parent,QString("select file"),dir,QString::fromStdString(files));
}
QString getchoosefolder(QWidget * parent,QString dir ){
    QFileDialog qfd(parent);
    qfd.setFileMode(QFileDialog::DirectoryOnly);
    qfd.setOption(QFileDialog::ShowDirsOnly,false);
    qfd.setWindowTitle("selct dir");
    if(!dir.isEmpty())qfd.setDirectory(dir);
    if(qfd.exec()==QDialog::Accepted){
        return qfd.selectedFiles().at(0);
    }
    return NULL;
}

QStringList getfilesfromfolder(string folder,string ext){
    QDir  dir(QString::fromStdString(folder));

    if(!dir.exists())
    {
        return QStringList();
    }
    QString qs=QString::fromStdString(ext);
    return dir.entryList(QStringList(qs));
    //return dir.entryList( );
}


int qgs_showimg( QGraphicsScene * qgs , QImage qi,int newimg,QGraphicsPixmapItem * oldqpgi){

    QGraphicsPixmapItem * pqgpi;
    if(newimg){
        pqgpi=new QGraphicsPixmapItem();
    }
    else{
        pqgpi=oldqpgi;
    }

    pqgpi->setPixmap(QPixmap::fromImage(qi));
    //pqgpi->setFlags(QGraphicsItem::GraphicsItemFlag::ItemIsMovable|QGraphicsItem::GraphicsItemFlag::ItemIsSelectable);//|QGraphicsItem::ItemEnabledChange|QGraphicsItem::ItemVisibleChange

    if(newimg){
        qgs->addItem(pqgpi);
    }

    //ui->graphicsView->setScene(&qgs);

    return 0;
}
#ifdef MY_USE_OPENCV
int qgs_showimg( QGraphicsScene * qgs, cv::Mat m,int c,int newimg,QGraphicsPixmapItem * oldqgpi){
    QImage qi;
    mat2qimage(m,qi,c);
    return qgs_showimg(qgs,qi,newimg,oldqgpi);

}
#endif
void qgv_fitsize(QGraphicsView * graphicsView,QGraphicsPixmapItem* qgpi){
    QImage qi=qgpi->pixmap().toImage();

    double s=1;
    //cout<<ui->graphicsView->height()<<endl;
    if((double)qi.height()/qi.width()>(double)graphicsView->height()/graphicsView->width()){
        s=1.0*(graphicsView->height())/qi.height();
    }
    else{
        s=1.0*(graphicsView->width())/qi.width();
    }

    graphicsView->setTransform(QTransform::fromScale(s,s));
}

#endif

#ifdef MY_USE_MATLAB
int Mat2mxArray(Mat const & imgSrc,mxArray * & pMat){
    int h,w,c;
    h = imgSrc.rows;
    w = imgSrc.cols;
    c = imgSrc.channels();
    int ttt=h*w*c;
    cout<<"converting Mat to mxArray size:"<<w<<"*"<<h<<"*"<<c<<" sum:"<<ttt<<"px"<<endl;

    // mxArray *pMat = NULL;
    unsigned char *input = NULL;
    if (c == 1) // gray image
    {
        mwSize dims[2] = {h, w};
        pMat = mxCreateNumericArray(2, dims, mxUINT8_CLASS, mxREAL);
        input = (unsigned char *)mxGetData(pMat);
        for (int i = 0; i < h; i++)
        {
            int tt;
            tt=0-h;
            for (int j = 0; j < w; j++)
            {	tt+=h;
                input[tt + i] =  (* imgSrc.row(i).col(j).data);
            }
        }
    }
    else if (c == 3) // 3-channel image
    {
        mwSize dims[3] = {h, w, c};

        pMat = mxCreateNumericArray(c, dims, mxUINT8_CLASS, mxREAL);

        input = (unsigned char *)mxGetData(pMat);
        int t=w*h;
        for (int i = 0; i < h; i++)
        {
            int tt=0-h;
            for (int j = 0; j < w; j++)
            {
                tt+=h;

                for (int ch = 0; ch < c; ch++)
                {
                    input[tt + i + t*ch] = (unsigned char)imgSrc.row(i).col(j).data[c - 1- ch];
                }
            }
        }
    }

    return 0;
}

int mxArray2Mat(const mxArray * mx,Mat & m){
    int dem=mxGetNumberOfDimensions(mx);
    int M=mxGetM(mx);
    int N=mxGetN(mx);
    int i,j;
    int t,tt;
    bool type=mxIsUint8(mx);
    if(type){
        if(dem>3){
            cout<<"too many demensions"<<endl;
            return -1;
        }
        uchar * imgdata=(uchar *)mxGetPr(mx);

        if(dem==3){
            N/=3;
            tt=M*N;
            m=Mat(M,N,CV_8UC3);
            for(i=0;i<M;i++){
                t=0-M;
                for(j=0;j<N;j++){
                    t+=M;
                    m.row(i).col(j).data[2]=imgdata[t+i];
                    m.row(i).col(j).data[1]=imgdata[t+i+tt];
                    m.row(i).col(j).data[0]=imgdata[t+i+tt+tt];
                }
            }

        }
        else{
            m=Mat(M,N,CV_8UC1);
            for(i=0;i<M;i++){
                t=0-M;
                for(j=0;j<N;j++){
                    t+=M;
                    m.row(i).col(j).data[0]=imgdata[t+i];
                }
            }



        }
    }
    else{
        cout<<"not uint8 can't convert"<<endl;
        return -1;

    }
    return 0;

}
#endif

#ifdef MY_USE_OPENCV
int mat2qimage(const cv::Mat& matt,QImage & qim,int c) {
    using namespace cv;
    cv::Mat mat=matt;
    if(mat.channels()==1){
        cv::Mat zero(mat.size(),mat.type(),cv::Scalar::all(0));

        std::vector<Mat> vm;
        if(c==0){
            vm.push_back(mat);
            vm.push_back(zero);
            vm.push_back(zero);
        }
        else if(c==1){
            vm.push_back(zero);vm.push_back(mat);
            vm.push_back(zero);

        }
        else if(c==2){
            vm.push_back(zero);
            vm.push_back(zero);vm.push_back(mat);

        }
        else{
            vm.push_back(mat);vm.push_back(mat);vm.push_back(mat);

        }
        cv::merge(vm,mat);

    }


    cv::Mat rgb;
    cvtColor(mat, rgb, CV_BGR2RGB);
    QImage tqi;

    tqi=QImage((const unsigned char*)(rgb.data), rgb.cols, rgb.rows,rgb.cols*3, QImage::Format_RGB888);
    qim=tqi.copy();
    return 0;
}

int qimage2mat(const QImage& qi,cv::Mat & matr) {
    using namespace cv ;
    cout<<(qi.format()==QImage::Format_RGB32)<<endl;
    QImage qimage=qi.copy();
    cv::Mat mat;
    if(qimage.format()==QImage::Format_RGB888){
        mat = cv::Mat(qimage.height(), qimage.width(), CV_8UC3, (uchar*)qimage.bits(), qimage.bytesPerLine());
    }
    else if(qimage.format()==QImage::Format_RGB32){
        mat = cv::Mat(qimage.height(), qimage.width(), CV_8UC4, (uchar*)qimage.bits(), qimage.bytesPerLine());
        const int from_to[] = { 0,0, 1,1, 2,2 };
        cv::Mat tm(mat.size(),CV_8UC3);
        mixChannels(&mat,1,&tm,1,from_to,3);
        mat=tm;
        //imshow("aaa",tm);
        cout<<tm.channels()<<endl;
    }
    else{
        cout<<(int)qimage.format()<<endl;
    }
    //cv::Mat mat2 = cv::Mat(mat.rows, mat.cols, CV_8UC3 );
    //int from_to[] = { 0,0, 1,1, 2,2 };
    //cv::mixChannels( &mat, 1, &mat2, 1, from_to, 3 );

    //Mat mat2=mat.clone();
    //cv::cvtColor(mat2,matr,CV_RGB2BGR);
    //imshow("t",matr);
    matr=mat;
    return 0;
}
#endif


#ifdef MY_USE_POINT_IO





//using namespace cv;
pointio::pointio()
{
}


pointset pointio::getpointsetformfile(string filename,int zindex)
{
    int xs=0;//sub points

    ifstream ifs;
    //filename="../3dbuild/"+filename;
    //cout<<"nowdir"<<QDir::currentPath().toStdString()<<"exist"<<QFile::exists(QString::fromStdString(filename))<<endl;
    ifs.open(filename.c_str());
    //cout<<"open="<<ifs.is_open();
    //cerr << "Error: " << strerror(errno)<<filename<<endl;
    point3d p,tp;
    tp.x=-1;
    char c;
    vector<point3d> pointlist;

    while(1){
        double t1,t2;
        //ifs>>p.x>>c>>p.y;
        ifs>>t1;
        ifs.get();
        //fgetchar();
        ifs>>t2;
        p.x=t1*expend;
        p.y=t2*expend;

        if(zindex){
            p.z=zindex;
        }
        else{
            ifs.get();
            double t3;
            ifs>>t3;
            p.z=t3;
        }

        p.z*=expend;

        if(ifs.fail()) break;


        if(processnose && tp.x!=-1){
            if(pointdistance( tp,p)>100||p.x<=0||p.y<=0){
            }
            else{
                if(xs++%((xsep+ysep)/2)==0)
                    pointlist.push_back(p);
            }
            tp=p;
        }
        else{
            pointlist.push_back(p);
        }




    }
    return pointlist;
}

void pointio::outputpoints(pointset pointlist)
{
    int i;
    int size=pointlist.size();
    for(i=0;i<size;i++){
        point3d p = pointlist.at(i);
        cout<<"x:"<<p.x<<" y:"<<p.y<<" z:"<<p.z<<endl;
    }
}

pointsetlist pointio::getpointsinfolder(string dirname,int autoz)
{
    int zs=0;
    QDir qd;
    pointsetlist psl;
    pointset ps;
    qd.setPath((QString)dirname.c_str());
    //qd.setCurrent();
    QStringList files= qd.entryList( QStringList(QString("*.txt")));
    int i;
    int len=files.size();
    for(i=0;i<len;i++){
        //cout<<files.at(i).toStdString()<<endl;
        string fn=dirname+QDir::separator().toLatin1()+files.at(i).toStdString();
        // cout<<fn<<endl;
        if(autoz){
            ps=getpointsetformfile(fn,autoz++);
        }
        else{
            ps=getpointsetformfile(fn);
        }
        if(zs++ %zsep==0)
            psl.push_back(ps);
    }
    //cout<<psl.size()<<endl;
    return psl;

}

int pointio::pointdistance(point3d p1, point3d p2)
{
    return fabs(p1.x-p2.x)+fabs(p1.y-p2.y)+fabs(p1.z-p2.z);
}
#ifdef MY_USE_OPENCV
void pointio::points2img(string dirnm,int image2d)
{
    using namespace cv;
    int fixx=0;
    int fixy=0;
    int imagew,imageh;
    imagew=imageh=525;
    // fixx=fixy=0;
    QDir qd;
    qd.setPath(QString(dirnm.c_str()));
    QStringList files=qd.entryList(QStringList(QString("*.txt")));
    int len=files.length(),lenj;
    int i,j;
    cv::Mat m;
    pointset ps;
    for(i=0;i<len;i++){
        string fn=dirnm+QDir::separator().toLatin1()+files.at(i).toStdString();
        ps=getpointsetformfile(fn,image2d);
        lenj=ps.size();
        m=Mat(imageh,imagew,CV_8UC3,Scalar::all(0));
        for(j=0;j<lenj;j++){
            mypoint p= ps.at(j);
            uchar * pi=m.ptr<uchar>(p.y+fixy,p.x+fixx);
            *(pi)=255;
            *(pi+1)=255;
            *(pi+2)=255;
        }
        imwrite(fn+".bmp",m);
        //break;

    }
}

#endif


#ifdef MY_USE_VTK


vtkSmartPointer<vtkPolyData> topolydata(vtkSmartPointer<vtkDataObject> src)
{
    vtkSmartPointer<vtkGeometryFilter> gf=vtkSmartPointer<vtkGeometryFilter>::New();
    gf->SetInputData(src);
    gf->Update();
    return  gf->GetOutput();
}

vtkSmartPointer<vtkPolyData> exactouterconn(vtkSmartPointer<vtkPolyData> src){

    vtkSmartPointer<vtkPolyDataConnectivityFilter> confilter=vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
    confilter->SetInputData(src);
    confilter->SetExtractionModeToLargestRegion();
    confilter->Update();
    vtkSmartPointer<vtkPolyData> lpd=vtkSmartPointer<vtkPolyData>::New();
    lpd->DeepCopy(confilter->GetOutput());

    vtkSmartPointer<vtkAppendPolyData> appd=vtkSmartPointer<vtkAppendPolyData>::New();
    appd->AddInputData(lpd);
    confilter->SetExtractionModeToAllRegions();
    confilter->Update();
    int rn=confilter->GetNumberOfExtractedRegions();

    confilter->SetExtractionModeToSpecifiedRegions();
    confilter->InitializeSpecifiedRegionList () ;

    vtkSmartPointer<vtkPolyData> tpl=vtkSmartPointer<vtkPolyData>::New();
    for(int i=0;i<rn;i++){
        confilter->AddSpecifiedRegion(i);
        if(i>0) confilter->DeleteSpecifiedRegion(i-1);
        confilter->Update();
        if(!boundisin(lpd->GetBounds(),confilter->GetOutput()->GetBounds())){
            tpl->DeepCopy( confilter->GetOutput());
            appd->AddInputData(tpl);
        }
    }
    //    confilter->Update();
    appd->Update();
    cout<<appd->GetOutput()->GetNumberOfPoints() <<"  " << appd->GetOutput()->GetNumberOfLines()<<endl;
    vtkSmartPointer<vtkPolyData> rst=vtkSmartPointer<vtkPolyData>::New();
    rst->DeepCopy(appd->GetOutput());


    return rst;
}

//check if sudully change lenght
static double __checkisinner_var_=-1;
static int checkisinner(double a,double b){
    return 1;
    double & pre=__checkisinner_var_;
    double len=abs(a-b);
    if(pre<0){
        pre=len;
        return 1;
    }
    //
    if(pre-len>1){
        // cout<<"drop"<<pre<<" "<<len<<endl;
        return 0;
    }
    else{

        // cout<<"keep"<<pre<<" "<<len<<endl;
        pre=len;
    }
    return 1;
}
//check if sudully change lenght
static double __checkisinner_var2_=-1;
static int checkisinner2(double a,double b){
    return 1;
    double & pre=__checkisinner_var_;
    double len=abs(a-b);
    if(pre<0){
        pre=len;
        return 1;
    }
    //
    if(pre-len>1){
        cout<<"drop"<<pre<<" "<<len<<endl;
        return 0;
    }
    else{

        // cout<<"keep"<<pre<<" "<<len<<endl;
        pre=len;
    }
    return 1;
}

static vtkSmartPointer<vtkPoints> combinne4(vtkSmartPointer<vtkPoints> p1,vtkSmartPointer<vtkPoints> p2,vtkSmartPointer<vtkPoints> p3,vtkSmartPointer<vtkPoints> p4,double distance=5){
    vtkSmartPointer<vtkKdTreePointLocator> ptrl1=vtkSmartPointer<vtkKdTreePointLocator>::New();
    vtkSmartPointer<vtkKdTreePointLocator> ptrl2=vtkSmartPointer<vtkKdTreePointLocator>::New();
    vtkSmartPointer<vtkKdTreePointLocator> ptrl3=vtkSmartPointer<vtkKdTreePointLocator>::New();
    vtkSmartPointer<vtkKdTreePointLocator> ptrl4=vtkSmartPointer<vtkKdTreePointLocator>::New();

    vtkSmartPointer<vtkPolyData> pl1=vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkPolyData> pl2=vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkPolyData> pl3=vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkPolyData> pl4=vtkSmartPointer<vtkPolyData>::New();

    pl1->SetPoints(p1);
    pl2->SetPoints(p2);
    pl3->SetPoints(p3);
    pl4->SetPoints(p4);


    ptrl1->SetDataSet((pl1));
    ptrl2->SetDataSet(pl2);
    ptrl3->SetDataSet(pl3);
    ptrl4->SetDataSet(pl4);

    ptrl1->BuildLocator();
    ptrl2->BuildLocator();
    ptrl3->BuildLocator();
    ptrl4->BuildLocator();


    vtkSmartPointer<vtkPoints> rst=vtkSmartPointer<vtkPoints>::New();



    vtkIdType i,pi;
    double *p,*nearp;
    double d2;
    vtkSmartPointer< vtkIdList> vids =vtkIdList::New();
    for(i=0;i<p1->GetNumberOfPoints();i++){
        p= p1->GetPoint(i);

        ptrl1->FindPointsWithinRadius(20,p,vids);
        if(vids->GetNumberOfIds()>45){
            rst->InsertNextPoint(p);
            continue;
        }

        pi=ptrl2->FindClosestPoint(p);
        nearp=p2->GetPoint(pi);
        d2=vtkMath::Distance2BetweenPoints(p,nearp);
        if(d2<distance){
            rst->InsertNextPoint(p);
            //            ptrl1->FindPointsWithinRadius(30,p,vids.GetPointer());
            //            cout<<"point nearby count"<<vids->GetNumberOfIds()<<endl;
            continue;
        }

        pi=ptrl3->FindClosestPoint(p);
        nearp=p3->GetPoint(pi);
        d2=vtkMath::Distance2BetweenPoints(p,nearp);
        if(d2<distance){
            rst->InsertNextPoint(p);
            continue;
        }

        pi=ptrl4->FindClosestPoint(p);
        nearp=p4->GetPoint(pi);
        d2=vtkMath::Distance2BetweenPoints(p,nearp);
        if(d2<distance){
            rst->InsertNextPoint(p);
            continue;
        }

    }


    for(i=0;i<p2->GetNumberOfPoints();i++){
        p= p2->GetPoint(i);


        ptrl2->FindPointsWithinRadius(30,p,vids);
        if(vids->GetNumberOfIds()>70){
            rst->InsertNextPoint(p);
            continue;
        }

        pi=ptrl3->FindClosestPoint(p);
        nearp=p3->GetPoint(pi);
        d2=vtkMath::Distance2BetweenPoints(p,nearp);
        if(d2<distance){
            rst->InsertNextPoint(p);
            continue;
        }

        pi=ptrl4->FindClosestPoint(p);
        nearp=p4->GetPoint(pi);
        d2=vtkMath::Distance2BetweenPoints(p,nearp);
        if(d2<distance){
            rst->InsertNextPoint(p);
            continue;
        }

    }


    for(i=0;i<p3->GetNumberOfPoints();i++){
        p= p3->GetPoint(i);

        ptrl3->FindPointsWithinRadius(30,p,vids);
        if(vids->GetNumberOfIds()>70){
            rst->InsertNextPoint(p);
            continue;
        }

        pi=ptrl4->FindClosestPoint(p);
        nearp=p4->GetPoint(pi);
        d2=vtkMath::Distance2BetweenPoints(p,nearp);
        if(d2<distance){
            rst->InsertNextPoint(p);
            continue;
        }

    }


    for(i=0;i<p4->GetNumberOfPoints();i++){
        p= p4->GetPoint(i);

        ptrl4->FindPointsWithinRadius(30,p,vids);
        if(vids->GetNumberOfIds()>70){
            rst->InsertNextPoint(p);
            continue;
        }

    }

    cout<<"after conbine "<<rst->GetNumberOfPoints()<<endl;
    return rst;
}

int pointisinbounds(double * p,double * bounds){
    double *b=bounds;
    return p[0]>=b[0] && p[0]<=b[1] && p[1]>=b[2] && p[1]<=b[3] && p[2]>=b[4] && p[2]<=b[5];
}

void processpointspieaces(vector<vtkSmartPointer<vtkPointSet> > &psl, int maxdis)
{
    if(psl.size()<1) return ;
    int maxlen=maxdis*maxdis;
    vector<vtkSmartPointer<vtkPointSet> > rst;
    vtkSmartPointer<vtkKdTreePointLocator> ptrl1=vtkSmartPointer<vtkKdTreePointLocator>::New();
    vtkSmartPointer<vtkKdTreePointLocator> ptrl2=vtkSmartPointer<vtkKdTreePointLocator>::New();
    vtkSmartPointer<vtkKdTreePointLocator> ptrl3=vtkSmartPointer<vtkKdTreePointLocator>::New();

    vtkSmartPointer< vtkPointSet> ps;//=psl.at(0);
    vtkSmartPointer< vtkPointSet> ps1;
    vtkSmartPointer<vtkPointSet> ps2;

    int i,j;
    double p[3];
    double length=0;
    vtkIdType vid,vid2;
    vtkSmartPointer<vtkPointSet> trst;
    for(i=0;i<psl.size();i++){
        trst=vtkSmartPointer<vtkPolyData>::New();
        trst->SetPoints(vtkSmartPointer<vtkPoints>::New());
        if(i==0){
            vid=-1;
        }
        else {
            ps=psl.at(i-1);
            ptrl1=vtkSmartPointer<vtkKdTreePointLocator>::New();
            ptrl1->SetDataSet(ps);

            ptrl1->BuildLocator();
        }

        if(i==psl.size()-1){
            vid=-1;
        }
        else{
            ps2=psl.at(i+1);
            ptrl3=vtkSmartPointer<vtkKdTreePointLocator>::New();
            ptrl3->SetDataSet(ps2);

            ptrl3->BuildLocator();
        }

        ps1=psl.at(i);
        for(j=0;j<ps1->GetNumberOfPoints();j++){
            ps1->GetPoint(j,p);
            if(i!=0)
                vid=   ptrl1->FindClosestPointWithinRadius(maxlen,p,length);
            if(i!=psl.size()-1)
                vid2=   ptrl3->FindClosestPointWithinRadius(maxlen,p,length);
            if(vid!=-1||vid2!=-1){
                trst->GetPoints()->InsertNextPoint(p);
            }
        }
        cout<<trst->GetNumberOfPoints()<<" : "<<ps1->GetNumberOfPoints()<<" : "<<i<<endl;
        rst.push_back(trst);
    }
    psl=rst;
}

void insertpointspieaces(vector<vtkSmartPointer<vtkPointSet> > &psl)
{
    if(psl.size()<1) return ;

    int insertlen=8;

    vector<vtkSmartPointer<vtkPointSet> > rst;
    vtkSmartPointer<vtkKdTreePointLocator> ptrl1=vtkSmartPointer<vtkKdTreePointLocator>::New();
    vtkSmartPointer<vtkKdTreePointLocator> ptrl2=vtkSmartPointer<vtkKdTreePointLocator>::New();
    vtkSmartPointer<vtkKdTreePointLocator> ptrl3=vtkSmartPointer<vtkKdTreePointLocator>::New();

    vtkSmartPointer< vtkPointSet> ps;//=psl.at(0);
    vtkSmartPointer< vtkPointSet> ps1;
    vtkSmartPointer<vtkPointSet> ps2;

    int i,j;
    double p[3],pp[3],tp[3];
    double length=0;
    vtkIdType vid,vid2;
    vtkSmartPointer<vtkPointSet> trst;
    int tcount=0,tcount2=0;
    for(i=0;i<psl.size();i++){
        tcount=tcount2=0;
        trst=vtkSmartPointer<vtkPolyData>::New();
        trst->SetPoints(vtkSmartPointer<vtkPoints>::New());
        if(i==0){
            vid=-1;
        }
        else {
            ps=psl.at(i-1);
            ptrl1=vtkSmartPointer<vtkKdTreePointLocator>::New();
            ptrl1->SetDataSet(ps);

            ptrl1->BuildLocator();
        }

        if(i==psl.size()-1){
            vid=-1;
        }
        else{
            ps2=psl.at(i+1);
            ptrl3=vtkSmartPointer<vtkKdTreePointLocator>::New();
            ptrl3->SetDataSet(ps2);

            ptrl3->BuildLocator();
        }

        ps1=psl.at(i);
        int count=ps1->GetNumberOfPoints();

        for(j=0;j<count;j++){
            ps1->GetPoint(j,p);
            if(i!=0){
                vid=ptrl1->FindClosestPoint(p);
                ps->GetPoint(vid,pp);
                double len=vtkMath::Distance2BetweenPoints(p,pp);
                if(len>insertlen){
                    tcount++;
                    int num= floor( sqrt(len-1));
                    tcount+=num-1;
//                     cout<<"insert point from :"<<endl;
//                     printarray(p,3);
//                     cout<<"to"<<endl;
//                     printarray(pp,3);
//                     cout<<endl;

                    tp[0]=p[0];tp[1]=p[1];tp[2]=p[2];
                    for (int i=1;i<num;i++){
                         tp[0]+=(pp[0]-p[0])/num;
                         tp[1]+=(pp[1]-p[1])/num;
                         tp[2]+=(pp[2]-p[2])/num;
                         ps1->GetPoints()->InsertNextPoint(tp);

                         //printarray(tp,3);
                    }

                }

            }

            if(i!=psl.size()-1){
                vid2=   ptrl3->FindClosestPoint(p);
                ps2->GetPoint(vid2,pp);
                double len=vtkMath::Distance2BetweenPoints(p,pp);
                if(len>insertlen){
                    int num= floor( sqrt(len));
                    tcount++;
                    tcount2+=num-1;
//                     cout<<"insert point from :"<<endl;
//                     printarray(p,3);
//                     cout<<"to"<<endl;
//                     printarray(pp,3);
//                     cout<<endl;

                    tp[0]=p[0];tp[1]=p[1];tp[2]=p[2];
                    for (int i=1;i<num;i++){
                         tp[0]+=(pp[0]-p[0])/num;
                         tp[1]+=(pp[1]-p[1])/num;
                         tp[2]+=(pp[2]-p[2])/num;
                         ps1->GetPoints()->InsertNextPoint(tp);

                         //printarray(tp,3);
                    }

                }
            }

        }
        //fix chee
        if(i==0 && false){
            double * c=ps1->GetCenter();
            double * b=ps1->GetBounds();
            ps1->GetPoints()->InsertNextPoint(c);
            ps1->GetPoints()->InsertNextPoint( (c[0]+b[0])/2,(c[1]+b[2])/2,i );
            ps1->GetPoints()->InsertNextPoint( (c[0]+b[0])/2,(c[1]+b[3])/2,i );
            ps1->GetPoints()->InsertNextPoint( (c[0]+b[1])/2,(c[1]+b[2])/2,i );
            ps1->GetPoints()->InsertNextPoint( (c[0]+b[1])/2,(c[1]+b[3])/2,i );

            ps1->GetPoints()->InsertNextPoint( c[0]/2 + (c[0]+b[0])/2/2 , c[1]/2+ (c[1]+b[2])/2/2,i );
            ps1->GetPoints()->InsertNextPoint( c[0]/2+(c[0]+b[0])/2/2,c[1]/2+(c[1]+b[3])/2/2,i );
            ps1->GetPoints()->InsertNextPoint( c[0]/2+  (c[0]+b[1])/2/2,c[1]/2+(c[1]+b[2])/2/2,i );
            ps1->GetPoints()->InsertNextPoint( c[0]/2+(c[0]+b[1])/2/2,c[1]/2+(c[1]+b[3])/2/2,i );
        }
 cout<<i<<":"<<tcount<<":"<<tcount2<<endl;
    }

}

vtkSmartPointer<vtkPolyData> fixpolydatabybound(vtkSmartPointer<vtkPolyData> src, vtkSmartPointer<vtkPolyData> bound,bool inverse)
{
    vtkSmartPointer<vtkKdTreePointLocator> ptrl=vtkSmartPointer<vtkKdTreePointLocator>::New();
    ptrl->SetDataSet(bound);
    ptrl->BuildLocator();

    double * bounds=bound->GetBounds();
    double * p,*p2;
    int pid;
    //vector<int> rmids;

    vtkSmartPointer< vtkIdTypeArray> ids=vtkSmartPointer<vtkIdTypeArray>::New();


    vtkSmartPointer<vtkSelectionNode> selectionNode =
            vtkSmartPointer<vtkSelectionNode>::New();
    selectionNode->SetFieldType(vtkSelectionNode::POINT);
    selectionNode->SetContentType(vtkSelectionNode::INDICES);

    for(int i=0;i<src->GetNumberOfPoints();i++){
        p=src->GetPoint(i);
        if(pointisinbounds(p,bounds)){
            pid= ptrl->FindClosestPoint(p);
            p2=bound->GetPoint(pid);
            double len=vtkMath::Distance2BetweenPoints(p,p2);
            if(len<5){
                //src->DeletePoint(pid);
                // i--;
                //  rmids.push_back(i);
                ids->InsertNextValue(i);
            }
        }
        else{
            ids->InsertNextValue(i);
        }
    }
    //    for(int i=0;i<rmids.size();i++){
    //        src->DeletePoint(rmids[i]);
    //    }
    //    src->RemoveDeletedCells();

    selectionNode->SetSelectionList(ids);
    selectionNode->GetProperties()->Set(vtkSelectionNode::CONTAINING_CELLS(), 1);

    if(inverse){
        // Get points that are NOT in the selection
        selectionNode->GetProperties()->Set(vtkSelectionNode::INVERSE(), 1); //invert the selection
    }
    vtkSmartPointer<vtkSelection> selection =
            vtkSmartPointer<vtkSelection>::New();
    selection->AddNode(selectionNode);

    vtkSmartPointer<vtkExtractSelection> extractSelection =
            vtkSmartPointer<vtkExtractSelection>::New();

    extractSelection->SetInputData(0, src);
#if VTK_MAJOR_VERSION <= 5
    extractSelection->SetInput(1, selection);
#else
    extractSelection->SetInputData(1, selection);
#endif
    extractSelection->Update();



    src=  topolydata(  extractSelection->GetOutput());
    return src;
}


vtkSmartPointer<vtkPolyData> exactsurface2d(vtkSmartPointer<vtkPolyData> src,   double stemp){
    bool enabledrop=false;

    double * bounds;// (xmin,xmax, ymin,ymax, zmin,zmax).
    double * center;
    int xmin=0,xmax=1,ymin=2,ymax=3,zmin=4,zmax=5;
    bounds=src->GetBounds();

    center=src->GetCenter();

    printarray(bounds,6);



    if(bounds[ymax]-bounds[ymin]<1e-300){
        vtkSmartPointer<vtkPoints> points=vtkSmartPointer<vtkPoints>::New();

        vtkSmartPointer<vtkPoints>po1=vtkSmartPointer<vtkPoints>::New();
        vtkSmartPointer<vtkPoints>po2=vtkSmartPointer<vtkPoints>::New();
        vtkSmartPointer<vtkPoints>po3 =vtkSmartPointer<vtkPoints>::New();
        vtkSmartPointer<vtkPoints>po4=vtkSmartPointer<vtkPoints>::New();

        vtkSmartPointer<vtkCellArray> ver=vtkSmartPointer<vtkCellArray>::New();
        int vernum=0;

        double i;
        vtkSmartPointer<vtkPlane> plan=vtkSmartPointer<vtkPlane>::New();
        vtkSmartPointer<vtkCutter> cutter=vtkSmartPointer<vtkCutter>::New();
        cutter->SetInputData(src);
        cutter->SetCutFunction(plan);
        plan->SetNormal(1,0,0);
        double * b2;
        for(i=bounds[xmin];i<=bounds[xmax];i+=stemp){
            plan->SetOrigin(i,bounds[ymin],bounds[zmin]);
            cutter->Update();

            if( cutter->GetOutput()->GetNumberOfPoints()>0){
                b2=cutter->GetOutput()->GetBounds();

                if(checkisinner(b2[xmin],center[0])){
                    points->InsertNextPoint(b2[xmin],b2[ymin],b2[zmin]);
                    ver->InsertNextCell(1,&vernum);
                    vernum++;

                }
                if(checkisinner(center[0],b2[xmax])){
                    points->InsertNextPoint(b2[xmax],b2[ymax],b2[zmax]);
                    ver->InsertNextCell(1,&vernum);
                    vernum++;
                }

                po1->DeepCopy(points);

            }
        }

        plan->SetNormal(0,0,1);
        __checkisinner_var_=__checkisinner_var2_=-1;
        for(i=bounds[zmin];i<=bounds[zmax];i+=stemp){
            plan->SetOrigin(bounds[xmin],bounds[ymin],i);
            cutter->Update();

            if( cutter->GetOutput()->GetNumberOfPoints()>0){

                b2=cutter->GetOutput()->GetBounds();
                if(checkisinner(b2[zmin],center[2])) {
                    points->InsertNextPoint(b2[xmin],b2[ymin],b2[zmin]);
                    po2->InsertNextPoint(b2[xmin],b2[ymin],b2[zmin]);
                    ver->InsertNextCell(1,&vernum);
                    vernum++;

                }
                if(checkisinner(b2[zmax],center[2])) {
                    points->InsertNextPoint(b2[xmax],b2[ymax],b2[zmax]);
                    po2->InsertNextPoint(b2[xmax],b2[ymax],b2[zmax]);
                    ver->InsertNextCell(1,&vernum);
                    vernum++;
                }


            }
        }

        plan->SetNormal(1,0,1);
        __checkisinner_var_=__checkisinner_var2_=-1;
        //int j;
        // double endp=max(bounds[xmax],bounds[zmax]);
        // double startp=min(bounds[xmin],bounds[zmin]);

        double lenx_y_z=bounds[zmax]-bounds[zmin]+(bounds[xmax]-bounds[xmin]);
        //normal chuizhi surface
        for(i=bounds[zmin];i<= bounds[zmin] + lenx_y_z ;i+=stemp){
            plan->SetOrigin(bounds[xmin],bounds[ymin],i);
            cutter->Update();

            if( cutter->GetOutput()->GetNumberOfPoints()>0){
                b2=cutter->GetOutput()->GetBounds();
                if(checkisinner(b2[zmin],center[2]) ){
                    points->InsertNextPoint(b2[xmin],b2[ymin],b2[zmax]);
                    po3->InsertNextPoint(b2[xmin],b2[ymin],b2[zmax]);
                    ver->InsertNextCell(1,&vernum);
                    vernum++;

                }
                if(checkisinner(center[2],b2[zmax])) {
                    points->InsertNextPoint(b2[xmax],b2[ymax],b2[zmin]);//fix
                    po3->InsertNextPoint(b2[xmax],b2[ymax],b2[zmin]);//fix
                    ver->InsertNextCell(1,&vernum);
                    vernum++;
                }

            }
        }


        plan->SetNormal(1,0,-1);
        __checkisinner_var_=-1;
        double * center=src->GetCenter();

        //         double tlen=max(bounds[xmax]-bounds[xmin],bounds[zmax]-bounds[zmin]);
        //         startp=center[0]-tlen/2;
        //        endp=center[0]+tlen/2;
        for(i=bounds[zmin];i<=bounds[zmin]+lenx_y_z;i+=stemp){
            plan->SetOrigin(bounds[xmax],bounds[ymin],i);
            cutter->Update();

            if( cutter->GetOutput()->GetNumberOfPoints()>0){
                b2=cutter->GetOutput()->GetBounds();
                if(checkisinner(center[2],b2[zmax])){
                    points->InsertNextPoint(b2[xmin],b2[ymin],b2[zmin]);
                    po4->InsertNextPoint(b2[xmin],b2[ymin],b2[zmin]);
                    ver->InsertNextCell(1,&vernum);
                    vernum++;

                }

                if(checkisinner(b2[zmin],center[2])){
                    points->InsertNextPoint(b2[xmax],b2[ymax],b2[zmax]);
                    po4->InsertNextPoint(b2[xmax],b2[ymax],b2[zmax]);
                    ver->InsertNextCell(1,&vernum);
                    vernum++;
                }

            }
        }
        //!!!!!!
        points= combinne4(po1,po2,po3,po4);
        vtkSmartPointer<vtkPolyData> tpoly=vtkSmartPointer<vtkPolyData>::New();
        tpoly=points2polywithver(points);
        vtkSmartPointer<vtkCleanPolyData> clean=vtkSmartPointer<vtkCleanPolyData>::New();

        //tpoly->SetPoints(points);
        //tpoly->SetVerts(ver);


        clean->SetInputData(tpoly);

        //        clean->set

        clean->Update();
        vtkSmartPointer<vtkPolyData> rst=vtkSmartPointer<vtkPolyData>::New();
        rst=clean->GetOutput();
        // rst=points2polywithver(rst->GetPoints());


        //
        //        savevtk(po1,"po1.vtk");
        //        savevtk(po2,"po2.vtk");
        //       savevtk(po3,"po3.vtk");
        // savevtk(clean->GetOutput(),"por.vtk");

        //another way
        if(true) return rst;


        cout<<tpoly->GetNumberOfPoints()<<" exact srucface rst:"<<clean->GetOutput()->GetNumberOfPoints()<<endl;



        //    vtkSmartPointer<vtkPolyDataConnectivityFilter> cf2=vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
        //    cf2->SetClosestPoint(0,0,0);

        vtkSmartPointer<vtkDelaunay3D> vd3d=vtkSmartPointer<vtkDelaunay3D>::New();

        vd3d->SetInputData(clean->GetOutput());
        vd3d->SetAlpha(5);
        //vd2d->SetSourceData();

        vd3d->Update();



        rst=topolydata( vd3d->GetOutput());

        //     confilter->SetInputData( );
        //     confilter->SetExtractionModeToAllRegions();
        //     confilter->Update();
        //     rn=confilter->GetNumberOfExtractedRegions();
        //    cout<<"rjionnum"<< rn<<endl;

        cout<<rst->GetNumberOfPoints()<<" final  rst"<<rst->GetNumberOfLines()<<endl;
        rst=exactouterconn(rst);

        cout<<rst->GetNumberOfPoints()<<" final  rst"<<rst->GetNumberOfLines()<<endl;



        return rst;



        //return clean->GetOutput();
        //      return tpoly;
    }
    else{
        assert(false);
    }


}

vtkSmartPointer<vtkPoints> pointio::mypointset2pointset(pointset pset){
    vtkSmartPointer<vtkPoints> points=vtkSmartPointer<vtkPoints>::New();
    int len=pset.size();
    for(int j=0;j<len;j++){
        point3d p=pset.at(j );
        points->InsertNextPoint(p.x,p.y,p.z);
    }
    return points;
}

void pointio::savevtkpointsets(vector<vtkSmartPointer<vtkPointSet> > pointsets)
{
    bool savez=false;
    QString fnm("output_ps");
    QDir qd(fnm);
    if(!qd.exists()) qd.mkpath(".");
    int i,j;
    char buffer[81];
    double *p;
    ofstream ofs;
    vtkSmartPointer< vtkPointSet> ps;
    for(i=0;i<pointsets.size();i++){
        ps=pointsets.at(i);
        sprintf(buffer,"%s%cp_%03d.txt",fnm.toStdString().c_str(),QDir::separator().toLatin1(),i);
        //cout<<buffer<<endl;

        ofs.open(buffer);

        for(j=0;j<ps->GetNumberOfPoints();j++){
            p=ps->GetPoints()-> GetPoint(j);
            ofs<<p[0]<<" "<<p[1];
            if(savez) ofs<<" "<<p[2];
            ofs<<endl;

        }
        ofs.close();
    }
}



vtkSmartPointer<vtkPolyData> readvtk(const char * filename){
    vtkSmartPointer <vtkGenericDataObjectReader> vtkgr=vtkSmartPointer<vtkGenericDataObjectReader>::New();
    vtkgr->SetFileName(filename);
    vtkgr->Update();
    vtkSmartPointer<vtkPolyData> poly=vtkSmartPointer<vtkPolyData>::New();
    if(vtkgr->IsFilePolyData()){
        poly=vtkgr->GetPolyDataOutput();
    }
    return poly;
}

vtkSmartPointer<vtkPolyData> removescalar(vtkSmartPointer<vtkPolyData> poly){
    poly->GetPointData()->SetScalars(0);
    poly->GetPointData()->SetVectors(0);
    poly->GetPointData()->SetNormals(0);
    return poly;
}


void savevtk(vtkSmartPointer<vtkPoints> points, const char *filename){
    vtkSmartPointer<vtkPolyData>poly=vtkSmartPointer<vtkPolyData>::New();
    poly->SetPoints(points);
    savevtk(poly,filename);
}


void savevtk( vtkSmartPointer<vtkPolyData> poly,const char * filename){

    cout<<"savefile "<<poly->GetPoints()->GetNumberOfPoints()<<endl;
    //    vtkSmartPointer<vtkPolyData> poly2=vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkPolyDataWriter> plyw=vtkSmartPointer<vtkPolyDataWriter>::New();


    // poly2->SetPoints(poly->GetPoints());
    plyw->SetInputData(poly);

    //plyw->SetFileTypeToASCII();

    plyw->SetFileName(filename);
    //plyw->Write();
    plyw->Update();


}
vector< vtkSmartPointer<vtkPointSet> > pointio::pointlist2poinsetvector(pointsetlist pointlist)
{
    vector< vtkSmartPointer<vtkPointSet> > rst;
    int len=pointlist.size();
    int i,j;
    for(i=0;i<len;i++){
        vtkSmartPointer<vtkPointSet> points=vtkSmartPointer<vtkPolyData>::New();
        points->SetPoints(vtkSmartPointer<vtkPoints>::New());
        pointset ps=pointlist.at(i);
        for(j=0;j<ps.size();j++){
            point3d p=ps.at(j);
            points->GetPoints()->  InsertNextPoint(p.x,p.y,p.z);
        }

        vtkSmartPointer<vtkCleanPolyData> clean=vtkSmartPointer<vtkCleanPolyData>::New();
        clean->SetInputData(  points2polywithver( points->GetPoints()));
        clean->Update();
        rst.push_back(clean->GetOutput());
    }
    return rst;
}
vtkSmartPointer<vtkPoints> pointio::pointlist2poinset(pointsetlist psl){
    vtkSmartPointer<vtkPoints> points=vtkSmartPointer<vtkPoints>::New();

    //points->InsertNextPoint();

    int leni=psl.size(),i,j;

    int lenj;
    for(i=0;i<leni;i++){
        lenj=psl.at(i).size();
        //cout<<"leni="<<lenj<<endl;
        for(j=0;j<lenj;j++){

            point3d p= psl.at(i).at(j);
            points->InsertNextPoint(p.x,p.y,p.z);
        }
    }
    return points;
}



vtkSmartPointer<vtkPolyData> points2polywithver(vtkSmartPointer<vtkPoints> points)
{
    vtkSmartPointer<vtkPolyData> poly=vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkCellArray> ver=vtkSmartPointer<vtkCellArray>::New();

    int i=0;
    for(i=0;i<points->GetNumberOfPoints();i++){
        ver->InsertNextCell(1,&i);
    }
    poly->SetVerts(ver);
    poly->SetPoints(points);

    return poly;
}



vtkSmartPointer<vtkUnstructuredGrid> getunstructgrid(const char * filenm)
{
    pointio pio;
    pointsetlist psl= pio.getpointsinfolder(filenm);

    vtkSmartPointer<vtkPoints>vtkpoints;
    vtkSmartPointer<vtkPolyVertex>polyvertex;
    vtkSmartPointer<vtkUnstructuredGrid> grid;
    vtkpoints=vtkPoints::New();
    polyvertex = vtkPolyVertex::New();
    grid=vtkUnstructuredGrid::New();

    int i,j,c=0;
    int leni=psl.size();

    int lenj;
    for(i=0;i<leni;i++){
        lenj=psl.at(i).size();
        //cout<<"leni="<<lenj<<endl;
        for(j=0;j<lenj;j++){

            point3d p= psl.at(i).at(j);
            vtkpoints->InsertPoint(c++,p.x,p.y,p.z+i*5);
        }
    }
    cout<<"points count"<<c<<endl;
    polyvertex->GetPointIds()->SetNumberOfIds(c);
    for(i=0;i<c;i++)//建立拓扑关系
    {
        polyvertex->GetPointIds()->SetId(i,i);
    }
    grid->SetPoints(vtkpoints);
    grid->InsertNextCell(polyvertex->GetCellType(),
                         polyvertex->GetPointIds());

    return grid;
    //map1->SetInputData(grid);
}

vtkSmartPointer<vtkPolyData> transform_back(vtkSmartPointer<vtkPoints> pt, vtkSmartPointer<vtkPolyData> pd);
vtkSmartPointer<vtkPolyData> getsurface(vtkSmartPointer<vtkPoints> points,int neiborsize,double space)
{
    vtkSmartPointer<vtkSurfaceReconstructionFilter> sfrf=vtkSurfaceReconstructionFilter::New();
    if(neiborsize>0)sfrf->SetNeighborhoodSize(neiborsize);//default is 20
    if(space>0)sfrf->SetSampleSpacing(space);
    //sfrf->SetNeighborhoodSize(0);
    //sfrf->SetSampleSpacing(0);

    //vtkSmartPointer<vtkUnstructuredGrid> vtkpd=vtkUnstructuredGrid::New();
    vtkSmartPointer<vtkPolyData> vtkpd=vtkPolyData::New();
    vtkpd->SetPoints( points);

    //is vpda needed?;
    //    vtkSmartPointer<vtkUnstructuredGridAlgorithm> vpda=vtkSmartPointer< vtkUnstructuredGridAlgorithm>::New();
    //    vpda->SetInputData(vtkpd);
    //    //cout<<vtkpd->GetPointData()
    //    vpda->Update();
    //    sfrf->SetInputConnection(vpda->GetOutputPort());
    sfrf->SetInputData(vtkpd);
    //sfrf->Update();

    vtkSmartPointer<vtkReverseSense> reverse=vtkReverseSense::New();
    vtkSmartPointer<vtkContourFilter> cf=vtkContourFilter::New();
    //cf->SetInputConnection();

    cf->SetInputConnection(sfrf->GetOutputPort());
    cf->SetValue(0,0.0);

    reverse->SetInputConnection(cf->GetOutputPort());
    reverse->ReverseCellsOn();
    reverse->ReverseNormalsOn();
    reverse->Update();

    vtkSmartPointer<vtkPolyData> newSurf = transform_back( points, reverse->GetOutput());
    return newSurf;

}

vtkSmartPointer< vtkPoints> projectPoints( vtkSmartPointer< vtkPoints> pt0,double direct[]){
    vtkSmartPointer< vtkPlane> plane = vtkSmartPointer< vtkPlane>::New();
    plane->SetOrigin(0,0,0);
    plane->SetNormal(direct);

    vtkSmartPointer< vtkPoints> newPts = vtkSmartPointer< vtkPoints>::New();
    double proj[3];
    for(int i=0;i<pt0->GetNumberOfPoints();i++){
        double p[3];
        pt0->GetPoint(i,p);
        plane->ProjectPoint(p, proj);

        newPts->InsertNextPoint(proj);
    }
    return newPts;
}


vtkSmartPointer<vtkPolyData> transform_back(vtkSmartPointer<vtkPoints> pt, vtkSmartPointer<vtkPolyData> pd)
{
    // The reconstructed surface is transformed back to where the
    // original points are. (Hopefully) it is only a similarity
    // transformation.

    // 1. Get bounding box of pt, get its minimum corner (left, bottom, least-z), at c0, pt_bounds

    // 2. Get bounding box of surface pd, get its minimum corner (left, bottom, least-z), at c1, pd_bounds

    // 3. compute scale as:
    //       scale = (pt_bounds[1] - pt_bounds[0])/(pd_bounds[1] - pd_bounds[0]);

    // 4. transform the surface by T := T(pt_bounds[0], [2], [4]).S(scale).T(-pd_bounds[0], -[2], -[4])



    // 1.
    double pt_bounds[6];  // (xmin,xmax, ymin,ymax, zmin,zmax)
    pt->GetBounds(pt_bounds);


    // 2.
    double pd_bounds[6];  // (xmin,xmax, ymin,ymax, zmin,zmax)
    pd->GetBounds(pd_bounds);

    //   // test, make sure it is isotropic
    //   std::cout<<(pt_bounds[1] - pt_bounds[0])/(pd_bounds[1] - pd_bounds[0])<<std::endl;
    //   std::cout<<(pt_bounds[3] - pt_bounds[2])/(pd_bounds[3] - pd_bounds[2])<<std::endl;
    //   std::cout<<(pt_bounds[5] - pt_bounds[4])/(pd_bounds[5] - pd_bounds[4])<<std::endl;
    //   // TEST


    // 3
    double scale = (pt_bounds[1] - pt_bounds[0])/(pd_bounds[1] - pd_bounds[0]);


    // 4.
    vtkSmartPointer<vtkTransform> transp = vtkSmartPointer<vtkTransform>::New();
    transp->Translate(pt_bounds[0], pt_bounds[2], pt_bounds[4]);
    transp->Scale(scale, scale, scale);
    transp->Translate(- pd_bounds[0], - pd_bounds[2], - pd_bounds[4]);

    vtkSmartPointer<vtkTransformPolyDataFilter> tpd = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
#if VTK_MAJOR_VERSION <= 5
    tpd->SetInput(pd);
#else
    tpd->SetInputData(pd);
#endif
    tpd->SetTransform(transp);
    tpd->Update();


    return tpd->GetOutput();
}

//#define MY_USE_PCL
#ifdef MY_USE_PCL

void getnomals(vtkSmartPointer<vtkPolyData> polydata,int redirect )
{

    //polydata->
    double * bounds=polydata->GetBounds();
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);


    vtkPolyData * const xp=polydata.GetPointer();
    //vtkPolyData* const xp=NULL;

    pcl::io::vtkPolyDataToPointCloud(xp,*cloud);
    //... read, pass in or create a point cloud ...

    // Create the normal estimation class, and pass the input dataset to it
    pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;
    ne.setInputCloud (cloud);

    // Create an empty kdtree representation, and pass it to the normal estimation object.
    // Its content will be filled inside the object, based on the given input dataset (as no other search surface is given).
    pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ());
    ne.setSearchMethod (tree);

    // Output datasets
    pcl::PointCloud<pcl::Normal>::Ptr cloud_normals (new pcl::PointCloud<pcl::Normal>);

    pcl::PointCloud<pcl::PointNormal>::Ptr ppn(new pcl::PointCloud<pcl::PointNormal>);




    // Use all neighbors in a sphere of radius 3cm
    //ne.setRadiusSearch (0.03);

    //ne.setRadiusSearch(3);
    //ne.setKSearch(12);
    ne.setKSearch(24);

    if(redirect)
        ne.setViewPoint((bounds[0]+bounds[1])/2,(bounds[1]+bounds[2])/2,(bounds[3]+bounds[4])/2);

    // Compute the features
    ne.compute (*cloud_normals);

    // cloud_normals->points.size () should have the same size as the input cloud->points.size ()*

    cout<<"got "<<cloud_normals->points.size()<<" nomals"<<endl;

    int i;
    pcl::PointNormal pn;
    pcl::PointXYZ p;
    pcl::Normal n;
    for(i=0;i<cloud_normals->points.size();i++){
        p=cloud->points.at(i);
        n=cloud_normals->points.at(i);
        //cout<<i<<" "<<n;
        pn.x=p.x;
        pn.y=p.y;
        pn.z=p.z;
        pn.normal_x=n.normal_x;
        pn.normal_y=n.normal_y;
        pn.normal_z=n.normal_z;
        ppn->push_back(pn);
    }


    pcl::io::pointCloudTovtkPolyData(*ppn,polydata.GetPointer());
}

#endif




#endif

#endif


#ifdef MY_USE_ITK

itk::MetaDataDictionary getdicomtags(string filename)
{

}

itk::Image<signed short,2>::Pointer readdicom(string filename)
{



    typedef signed short pixelType;
    const unsigned int dimension=2;
    typedef itk::Image<pixelType,dimension> imageType;

    typedef itk::ImageFileReader<imageType> readerType;

    //readerType::Pointer reader=readerType::New();
    itk::ImageFileReader<imageType>::Pointer reader=readerType::New();

    reader->SetFileName(filename);


    typedef itk::GDCMImageIO imageIoType;
    imageIoType::Pointer gdcmImageIo = imageIoType::New();
    reader->SetImageIO(gdcmImageIo);


    try{
        reader->Update();
    }
    catch(itk::ExceptionObject & e){
        std::cerr << "exception in file reader " << std::endl;
        std::cerr << e.GetDescription() << std::endl;
        std::cerr << e.GetLocation() << std::endl;

    }

    //imageType::Pointer image=reader->GetOutput();
    itk::Image<pixelType,dimension>::Pointer image=reader->GetOutput();
    return image;

}

void getdicomdata(itk::Image<signed short,2>::Pointer img,ostream & out){
    string pcolumns=searchtagfromimage(img,"","Columns");
    string prows=searchtagfromimage(img,"","Rows");
    int col=atoi(pcolumns.c_str());
    int row=atoi(prows.c_str());
    int i,j;
    out<<"img row="<<row<<"col="<<col<<endl;
    for(i=0;i<row;i++){
        for(j=0;j<col;j++){
            itk::Image<signed short,2>::IndexType pi={{i,j}};
            out<<img->GetPixel(pi)<<" ";
        }
        out<<endl;
    }

}


string searchtagfromimage(itk::Image<signed short,2>::Pointer image, string key, string label)
{
    if(key==""&&label=="") return "";
    typedef itk::MetaDataDictionary dictionaryType;

    dictionaryType & dictionary =image->GetMetaDataDictionary();

    typedef itk::MetaDataObject<std::string> metaDataStringType;



    itk::MetaDataDictionary::ConstIterator it= dictionary.Begin();
    while(it!=dictionary.End()){
        itk::MetaDataObjectBase::Pointer  entry = it->second;
        metaDataStringType::Pointer entryvalue =
                dynamic_cast<metaDataStringType *>( entry.GetPointer() );

        if( entryvalue ){
            std::string tagkey   = it->first;
            std::string tagvalue = entryvalue->GetMetaDataObjectValue();
            if(key!="" && tagkey==key) return tagvalue;
            std::string labelId;
            bool found =  itk::GDCMImageIO::GetLabelFromTag( tagkey, labelId );
            if( found )
            {
                if(labelId==label) return tagvalue;

            }


        }


        ++it;

    }

    return "";

}


void listtagsfromimage(itk::Image<signed short,2>::Pointer image,bool onlyimportant,ostream &out)
{

    typedef itk::MetaDataDictionary dictionaryType;

    dictionaryType & dictionary =image->GetMetaDataDictionary();

    typedef itk::MetaDataObject<std::string> metaDataStringType;



    itk::MetaDataDictionary::ConstIterator it= dictionary.Begin();
    while(it!=dictionary.End()){
        itk::MetaDataObjectBase::Pointer  entry = it->second;
        metaDataStringType::Pointer entryvalue =
                dynamic_cast<metaDataStringType *>( entry.GetPointer() );

        if( entryvalue ){
            std::string tagkey   = it->first;
            std::string tagvalue = entryvalue->GetMetaDataObjectValue();

            std::string labelId;
            bool found =  itk::GDCMImageIO::GetLabelFromTag( tagkey, labelId );
            if( found )
            {
                if(onlyimportant){

                    if( labelId=="Protocol Name" || labelId=="Image Position (Patient)"
                            || labelId=="Image Orientation (Patient)" || labelId=="Pixel Spacing"
                            ||labelId=="Slice Thickness" || labelId=="Rows"|| labelId=="Columns"
                            ||labelId=="Bits Allocated" || labelId=="Series Description"
                            ||labelId=="Patient's Name"||labelId=="Scanning Sequence" ||labelId=="Sequence Variant"
                            ||labelId=="Scan Options" || labelId=="Pixel Bandwidth" || labelId=="Acquisition Matrix"
                            ||labelId=="In-plane Phase Encoding Direction" ||labelId=="SAR" ||labelId=="dB/dt" ||labelId=="Patient Position"
                            ||labelId=="Series Number"||labelId=="Instance Number"||labelId=="Acquisition Number"
                            ||labelId=="Bits Stored"||labelId=="Window Center" || labelId=="Window Width"
                            ||labelId=="Window Center & Width Explanation" || labelId=="Study Description"
                            ||labelId=="Study ID"
                            ){
                        out << "(" << tagkey << ") " << labelId;
                        out << " = " << tagvalue.c_str() << std::endl;
                    }
                }
                else{
                    out << "(" << tagkey << ") " << labelId;
                    out << " = " << tagvalue.c_str() << std::endl;
                }
            }

            else if(!onlyimportant)
            {
                out << "(" << tagkey <<  ") " << "Unknown";
                out << " = " << tagvalue.c_str() << std::endl;
            }
        }


        ++it;

    }

}



void sepratedicom(QString sourcedir, QString destdir, string filter)
{
    QString dir=sourcedir;
    if(!dir.endsWith("/")) dir.append("/");
    QString ddir=destdir;
    QStringList qsl= getfilesfromfolder(dir.toStdString(),filter);
    int i=0;
    for(i=0;i<qsl.size();i++){
        QString fn=qsl.at(i);
        if(fn=="."||fn=="..")
            continue;
        fn=dir+fn;
        itk::Image<signed short,2>::Pointer image=readdicom(fn.toStdString());
        string pname=searchtagfromimage(image,"","Patient's Name");
        string pstudy=searchtagfromimage(image,"","Study Description");
        string pseries=searchtagfromimage(image,"","Series Description");
        string pseriesnumber= searchtagfromimage(image,"","Series Number");
        string pinstancenumber= searchtagfromimage(image,"","Instance Number");
        pinstancenumber.erase(pinstancenumber.find_last_not_of(" \n\r\t")+1);
        while (pinstancenumber.length()<5){
            pinstancenumber="0"+pinstancenumber;
        }

        string pprotocol= searchtagfromimage(image,"","Protocol Name");
        string newdir=ddir.toStdString()+pname+"/"+pstudy+"/"+pseries+"("+pseriesnumber+")/"+pprotocol+"/";
        QString qnewdir=QString::fromStdString(newdir);

        QDir qd(qnewdir);

        if(!qd.exists())
            cout<<qd.mkpath(qnewdir);



        QString qnewfilename=qnewdir+QString::fromStdString(pinstancenumber+"-")/*+ qsl.at(i)*/+QString(".dicom");
        QFile qf(fn);
        qf.copy(qnewfilename);
        //cout<<qf.errorString()<<endl;
        cout<<i<<":"<<fn<<" "<<qnewfilename<<endl;
    }
}

#include<itkDiscreteGaussianImageFilter.h>

itk::Image< signed short , 3 >::Pointer smoothitkimage(itk::Image< signed short , 3 >::Pointer img ){
    itk::DiscreteGaussianImageFilter< itk::Image< signed short , 3 >, itk::Image< signed short , 3 > >::Pointer filter=
            itk::DiscreteGaussianImageFilter< itk::Image< signed short , 3 >, itk::Image< signed short , 3 > >::New();
    filter->SetInput(img);
    filter->SetVariance(2);

    filter->Update();
    return filter->GetOutput();
}

#ifdef MY_USE_VTK

extern vtkSmartPointer<vtkImageData> readdicomfromfolder(string folder){
    itk::ImageSeriesReader< itk::Image< signed short , 3 > >::Pointer   readers=  itk::ImageSeriesReader< itk::Image< signed short , 3 > >::New() ;
    itk::GDCMImageIO::Pointer gdcmImageIO=itk::GDCMImageIO::New();
    readers->SetImageIO( gdcmImageIO );

    itk::GDCMSeriesFileNames::Pointer nameGenerator = itk::GDCMSeriesFileNames::New();

    nameGenerator->SetRecursive(true);
    //nameGenerator->SetUseSeriesDetails( true );
    // nameGenerator->AddSeriesRestriction("0008|0021" );

    nameGenerator->SetDirectory(folder);


    typedef std::vector< std::string >    SeriesIdContainer;

    const SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs();

    SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();
    SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();
    while( seriesItr != seriesEnd )
    {
        std::cout << seriesItr->c_str() << std::endl;
        ++seriesItr;
    }

    string seriesid=seriesUID.begin()->c_str();

    //typedef std::vector< std::string >   FileNamesContainer;
    itk::FilenamesContainer fileNames;

    fileNames=nameGenerator->GetFileNames(seriesid);

    cout<<"file count:"<<fileNames.size()<<endl;
    readers->SetFileNames(fileNames);

    readers->Update();
    //readers->UpdateLargestPossibleRegion();

    itk::ImageToVTKImageFilter<itk::Image<signed short,3> >::Pointer connector3d=itk::ImageToVTKImageFilter<itk::Image<signed short,3> >::New();

    //67 this shoud be a option
    if(1){
        itk::Image< signed short , 3 >::Pointer iimg =smoothitkimage( readers->GetOutput());
        connector3d->SetInput(iimg );
    }
    else{
        connector3d->SetInput(readers->GetOutput());
    }

    connector3d->Update();
    vtkSmartPointer<vtkImageData> vimgd=vtkSmartPointer<vtkImageData>::New();
    // vimgd=connector3d->GetOutput();
    vimgd->DeepCopy(connector3d->GetOutput());
    return vimgd;
}

#endif



#endif




